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Abstract 

We propose a generalization of the classical M/M/1 queue process. The resulting model is derived by applying 
fractional derivative operators to a system of difference-differential equations. This generalization includes both 
non-Markovian and Markovian properties, which naturally provide greater flexibility in modeling real queue 
systems than its classical counterpart. Algorithms to simulate M/M/1 queue process and the related linear 
birth-death process are provided. Closed-form expressions of the point and interval estimators of the parameters 
of these fractional stochastic models are also presented. These methods are necessary to make these models 
usable in practice. The proposed fractional M/M/1 queue model and the statistical methods are illustrated using 
S&P data. 

Keywords: Transient analysis. Fractional M/M/1 queue, Mittag-Leffler function. Fractional birth-death process. 
Parameter estimation. Simulation. 

1 Introduction 

The M/M/1 queue is without a doubt the simplest model for a queue length. It is characterized by arrivals determined 
by a Poisson process and an independent service time which is negative-exponentially distributed. It is relatively 
simple and yet the analysis of its transient behavior leads to considerable difficulties. The main source of these 
difficulties is the presence of a non- absorbing boundary at zero (empty queue). This means that the analysis becomes 
simpler when we consider models with absorbing boundaries. As a direct result, the state probability of a linear 
birth-death process, that is the probability that the queue length is n at a specific time has a particularly nice 
form. 

The aim of this paper is to study some related point processes governed by difference-differential equations containing 
fractional derivative operators. These processes are direct generalizations of the classical M/M/1 queue and the linear 
birth-death processes. It is well-known that a fractional derivative operator induces a non-Markovian behavior into a 
system [see 19]. Moreover, model parameter estimation and path generation algorithms of these fractional stochastic 
models are derived. Note that the proposed fractional point models (with Markovian and non-Markovian properties) 
are parsimonious, which is a desirable property for modeling real-world non-Markovian queueing systems. Also, 
observe that fractional point processes driven by fractional difference-differential equations such as the fractional 
Poisson, the fractional birth, the fractional death, and the fractional birth-death processes have been gaining 
attention more recently [see, e.g., 11, 4, 7, 14, 9]. 

The article is structured as follows. Section 2 presents the explicit construction of the fractional M/M/1 queue 
starting from the governing equations and a particular subordination relation. The main result derived in this section 
is the explicit form of the transient state probabilities for each value of the parameter of fractionality. Information 
regarding the steady-state behavior (stationary behavior) is also highlighted; in particular the fractional process 
shares the same steady-state behavior as the classical non-fractional case. In Section 3 we develop closed-form 
estimators (point and interval) for the model parameters in the case of the fractional linear birth-death process. 



This is preparatory for the similar subsequent analysis applied to the M/M/1 queue (Section 4). The article ends 
with an application which showed that our constructed estimators perform well in a real-world example. 



2 Results for a fractional process related to M/M/1 queues 

The classical M/M/1 queue process N{t), t > 0, that is the queue length in time can be described by the following 
difference-differential equations governing the state probabilities Pk{t) = Pr{N{t) = fc|A^(0) = i}, k > 0: 

JiPkit) = -(A + ^i)pk{t) + \pk-i{t) + tWk+i{t), fc > 1, 

f^Po{t) = -Xp„{t)+fipiit), (2.1) 
,Pfc(0) = Sk,i, 

where i G N U {0} is the initial number of individuals in the queue and S^^i is the Kronecker's delta. In (2.1) A > 
and /i > are the entrance and the service rates, respectively. 

To arrive at a possible fractional model we consider the Caputo fractional derivative D", a G (0, 1], with respect 
to time t. If p^{t) = Pr{7V"(i) = fc}, fc > 0, where N°'{t), is the fractional M/M/1 queue with parameter a, the 
generalized difference-differential equations for the state probabilities with arrival rate A > 0, service rate fi > and 
z > initial customers, read 

D?pnt) - -(A + M)Pfe(t) + Ap^l(t) + /iPfc+i(i), k > 1, 

D»p^it) = -Xp'^{t)+tipnt), (2.2) 

First, we will follow Bailey [2, 3] for the derivation of the probabilities p'^{t), k > 0, t > but adapting the method 
to take into considerations the presence of the Caputo derivative. The result obtained by Bailey is the so-called 
classical solution in terms of modified Bessel functions of the first kind. Note however that the derivation of the state 
probabilities in the classical case a = 1 can be carried out in several equivalent ways (see for example Champernowne 
[8], Parthasarathy [15], Abate and Whitt [1]). In the following we will first treat the solution derived by Bailey and 
then we will use a simpler but lesser known form due to Sharma [17]. 
We indicate G°'{z,t) — J^kLo '^''PkW ^ ^^"^ probability generating function. 

Theorem 2.1. The Laplace transform G°'{z^s) — e^**G'"(z, i) di, a e (0,1], can he written as 

^a, N - (1 - ^) [a2(g)]'+^ [1 - a2(g)]~^ I I ^ 1 »^ W 

G z,s =s — --^ —- , z <l, 3f? s >0. (2.3) 

where ai{s) and 02(5) are the zeros of f{z,s) — zs" — (1 — z)(^ — Xz). 
Proof. From (2.2), we can write 

[G"(z, t) - p^{t)] = -(A + A*) [G"(z, t) - f^{t)] + XzG^iz, t). (2.4) 
Using the equation on jjq (t) we have 

A"G"(z, t) = -AG"(z, t) - ^G"(z, t) + ^ip'^{t) + AzG"(z, t) + ^ [G"(z, t) ~ p'^{t)] , (2.5) 



and after simplifying, we obtain, for \z\ < 1, the Cauchy problem 

{zDfG^iz, i) = (1 - z) [G"(z, t){ii - Xz) - pLp^it)] 
lG"(z,0) = z\ 



(2.6) 
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Applying the Laplace transform G"{z, s) — e **G"(z, t) dt to (2.6) leads to 

z [s"G"(z, s) - s"-iG"(z, 0)1 (1 - z) [g"(z, s)(m - Az) - mPo (s) 
where p^(s) = e-'*p^{t) dt. After some simple algebraic calculations we then have 



G"(z,s) 



zs" - (1 - 2;)(/i - Az) 



\z\ < 1, 3fi(s) > 0. 



(2.7) 



(2.8) 



As the above function converges in \z\ < 1, the zeros of the numerator and the denominator should coincide. Let us 
indicate the zeros of the numerator as 

.s" + A + ± [(s" + A + Ai)2 - 4Am] 



ai2(s) = 

with \a2{s)\ < \ai{s)\, 3?(s) > 0. Note that 



2A 



(2.9) 



'ai(s) + 02(5) = (s" + A + /i)/A, 

ai(s)a2(s) = m/A, (2.10) 
[-A[l-a2(s)][l-ai(s)] =s°. 

By Rouche theorem [10, Page 168] we have that the only zero in the unit circle is 02(5). Therefore it follows that 



[a2(s)]'+^-M[l-a2(s)K(s)-0, 



which gives 



Po(s) 



[«2(.)] 



J+1 



/i[l-a2(s)] 
Now, by considering that 

zs" - (1 - z){fi - Xz) = ~X[z - ai(,s)] [z - 02(3)] , 
equation (2.8) can be rewritten as 

.a-l^^+'-(l-^) [a2is)Y^' [1-«2(.S)]-' 



G°(z,s) = s° 



-A [z - ai{s)] [z - a2{s)] 



(2.11) 
(2.12) 
(2.13) 

(2.14) 

□ 



In the following Theorem 2.2 we prove a subordination relation for the fractional queue N°'{t), i > 0, a G (0, 1]. 
This is essential for our next results. Before that, let us introduce some facts on the a-stable subordinator and its 
inverse process. 

Let us call V^^t), t>0, the a-stable subordinator (see for details Bertoin [5], cap. Ill) and let us define its inverse 
process as its hitting time 

E''{t)^M{s>0:V°'{s)>t}. (2.15) 
The processes V" {t) and i?" (i) are characterized by their Laplace transforms. For the a-stable subordinator we have 

Ee-«^° ^ e"*«°, ae(0,l], (2.16) 
and for its inverse process the time-Laplace transform reads 



/>oo 

/ e-«* (Pr{£;"(t) e ds}/ds) dt = C^e'"^" a e (0, 1]. 
Jo 



(2.17) 
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Theorem 2.2. Let N'^{t) = N{t), t > 0, be the classical M/M/l queue and let E°'{t), t > Q, a £ {0, 1], be an 
inverse a-stable subordinator (2.15) independent of N^{t). The fractional M/M/l queue N°'{t), t > 0, a G (0, 1], 
can be represented as 



N^'^t) = N^E^it)), t>0, ae(0,l), 
where the equality holds for the one- dimensional distribution. 
Proof. Let us consider the initial value problem 

izOfC^iz, i) = (1 - z) [G"(z, t){^ - \z) ~ pLp^it)] , 
[G"(z,0) = z\ 

which is equivalent to (2.2). Applying the Laplace transform we obtain 

z [s"G"(z,s)-s"-iG"(z, 0)1 = (1-z) \g"{z,s){ii-\z)- iip'^is) 
Note that if (2.18) holds we can write 



G"(z,s) 



■oo 



Pr{iV"(y) = k}¥v{E°'{t) E dy} 



dt 



POO POO 

/ e-^* / G(z,2/)Pr{i?"(t)edy} 
Jo Uo 



dt 



G{z,y). 



Q— 1 —ys'^ 



dy, 



and 



e'^'p^mt 



Po(y)Pr{i?"(i)Gdy} 



dt 



Po(y)s' 



a~l -ys" 



dy. 



We now show that (2.21) and (2.22) satisfy (2.20). Observe that 



G{z,y)e'y'"dyp-z' 



(l-z) 



(m-Az)/ Giz,y)e~y^''dy-^l p^{v)e-y'^ dy 



Consider the right hand side of (2.21). We can write 



(l-z) 



(m-Az) / G{z,y)e-y' dy - fi Poiy)e'y' dy 
Jo Jo 

e-y^" (1 - z) [{p - Xz)G{z, y) - pip^{y)] dy. 



Considering that G{z,y) and po{y) satisfy 

d 



z—-G{z, y) = (1 - z) [(^ - Az)G(z, y) - npo{y)] 
dy 



(2.18) 



(2.19) 



(2.20) 



(2.21) 



(2.22) 



(2.23) 



(2.24) 



(2.25) 
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we immediately obtain that 



(l-z) 



(m - A^) / G(z, y)e-y' dy - fi j poiy)e-y' dy 



(2.26) 



dy 



G{z,y)dy 



G{z,y)e 



y=o 

J 

Giz,y)e'y'"dy-z' 



G{z,y)e-y'"dy 



This concludes the proof. 



□ 



Using the Laplace transform (2.3) and the calculations carried out in Bailey [3] we can gain some insights on the 
mean value of the process. 



Theorem 2.3. We have that 



where 



T[a + 1) 



iy-tr-'fiy)dy, t>0, 



r(a) Jo 

is the Riemann-Liouville fractional integral [16]. 

Proof. By means of the Laplace transform (2.3) of the probability generating function we can write 

EiV"(s) = -^G"(z,s) 
dz 

[-X{z - ai){z - aa)] [{i + l)z' + al+\l - az)''] + A(2 ~ ai - aa) 
A2(l-ai)2(l-a2)2 
i + 1 + a2+^(l - 02)"^ 2-ai-a2 
A(l-ai)(l-a2) + A(l-ai)2(l-a2)2 

i + l + a^+^(l-a2)~^ 2A-(s" + A + m) 

s" ^ A2(l-ai)2(l-a2)2_ 

i + l + a^+^(l-a2"^) 2A-(s" + A + ^) 
1 + i + al+\l - a2)-^ X-^i s" 



(2.27) 



(2.28) 



(2.29) 



j_ ^ A - _^ a2+^(l - aa)^ 



X-fi s"-ia*2+'(l - 02) 



i , A-A* , Pq{s) 



s s 



a+l 
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Note the in the above calculation we have used the relation (2.10). Result (2.29) immediately yields 

E7V"(i) = z + (A - f,) ^° + /iJ>o" W- (2.30) 
1 (a + 1) 

□ 

Remark 2.1. The validity of Theorem 2.2 can he checked with the aid of formula (2.27) as follows. 

POO 

EN"{t)= EN\w)Pi{E"{t) G dw} (2.31) 
Jo 

= i + (A-^) / u.Pr{S"(<) e dw} + /i / pl{y)dyPr{E''{t) £dw}. 
Jo Jo Jo 

Therefore the time-Laplace transform, recalling that e~**dt Pr{£'"(t) G dw} — s°'~^e~^^ dw, can be written as 



EN"{s)^-+ ws"-ie-"'"°dw + /i / / pl{y)dys°"^e'""''' dw 
s Jo Jq Jo 



i A — u 

- + - 

s S' 



S S' 



a+1 



-ys 



/i- 



Now, by noticing that Pq{s°') / s —Pq{s)/s°' (see formula (2.12)j we arrive at 



s s"^^ s 



A Ao' 



2 



we arrive at 



i A — /i /iflj — Aai 



4+1 



(2.32) 



EiV«(s) = - + ^+M^, (2.33) 



which leads to formula (2.27). 

Remark 2.2. yl different form of formula (2.29) can be achieved by writing 

EN'^is) - - + ^ + ' ,^ °^ , (2.34) 
^ ' s s"+i s"(l - aa) 

_ i A-/i s"-iA4+\l - ai) 

where we used the fact that 1/(1 — 02) = — A(l — ai)/s". Furthermore, after considering 

s" + A + /i 

fli = > 02 = T — , (2.35) 
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Let us now address the problem of finding explicit results for the state probabilities = Pr{A^"(t) = fc|A^"(0) = i 
of the proposed fractional queue model. We start by using the subordination relation stated in Theorem (2.2) with 
the classical solution of the M/M/1 queue in terms of modified Bessel functions of the first kind. In the non-fractional 
case (a = 1) we have [3, Page 154] 



pUt) 



i(fe-i) 



/2, 



(2.37) 



(^+^)^ {a/,+,+2 (2(Am)^/V) - 2ix^^y/'h+k+l {2ix^iy/^l 

+Ai/.+fc(2(AAi)i/V)}dr, 



where /i/(z) is the modified Bessel function of the first kind. 

The state probabilities p'^{t), t > 0, k > 0, a G (0, 1] can thus be determined formally by subordination in the 
following way: 



/"OO 

Jo 

Using the time-Laplace transform p'^{s) — e^^'Vfe (^)di we have 

POO 

Pl{s)^ / p^(y)s"-ie-^^°dy 
Jo 



e-(A+M)yj^_^ {2{\^^Y/^y) s^-^e'y'"dy 

e-('+'''" {A/.+fc+2 (2(A/.)i/V) - 2(AM)^/'/.+fc+i (2(A/i)^/V) 
drj s'^-^e-y"^ dy 
-y(s-+x+,)j^_^ (2(Am)^/^2/) 

'e-(A+M)r (2(Am)'/V) - 2{\^Jif/^h+k+l (2(A/i)i 



+^Jih+k (2(A/i)i/V)} 

X\ 5(fe-2) 



S 



/2. 



(2(Am)^/V)}. 



dr / e-y"" dy, 



Applying the well-known Laplace transform for I,^{z) we get 



, ^ i(/c-2) 



2(Am) 



1/2 



k — i 



" + X + fi- J{s^ + X + fiy- [2(AAi)i/2l 



(s" + A + Ai)'- [2(AAi)i/2 



s 

OO 



A e-(«°+^+^)-W2 (2(Am)^/V) dr 
(Am)1/2 /"e-(^°+^+'^)-Wi (2(Am)^/V) 



(2.38) 



(2.39) 



(2.40) 
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/•OO 














-^^Ar 




S - 


X [(s° + A 





s" + A + /i - v/(s" + A + - 4A// 



-4A,]-/^ 
AV*^"'' 2(AAi)i/2 



2(A/i)i/2 

(j+fc+2) 



[(s" + A + - 4AAi; 



-1/2 



" + A + - + A + - 4A;u 



i+/£+2 



1/2 



-(i+fc+1) 



3" + A + ^ - + a + /^)2 - 4A/X 



i+fc + l 



X [(.s° + A + - 4A/i; 



-1/2 



i(fe-») 



2(Am) 



1/2 



-(i+fc) 



+ A + M - + A + /i)2 - 4A^ 



i+fe 



'1/2 



[(s" + A + A*)' " 4A/i] 

V(s« + A + Ai)2 - 4AAi 

(A/^)5('=-i)A 



-I i — k 



s"^(s" + A + - 4A/i 

(A/m)^(^-^)2(Aa.)V2 
s"V(s" + A + - 4A/i 

(A/M)^'^-^)/i 
+ A + - 4A/i 



s" + A + /i ~ ^(s" + A + /i)2 - 4A/i 
2(Am)i/2 

s° + A + /i - ^(s" + A + /i)2 - 4A^ 
2(Am)i/2 

s" + A + /i - ^(s" + A + /x)2 - 4A/I 
2(Am)i/2 

+ A + - ^(s" + A + ^)2 -4A^ 
2(A/i)i/2 



'i+fc+2 



i+fc+l 



Although the obtained Laplace transform p^(s) has a clear structure it cannot be inverted in a simple manner. Note 
anyway that it should be related to the Laplace transform of some generalizations of Bessel functions. 
In order to obtain more explicit results we must abandon the classical form of the state probabilities in terms of 
Bessel functions. We exploit instead a lesser known but certainly more appealing result due to Sharma [17, Chapter 
2]. In particular we refer to equation (2.2.16) at page 17 which we recall here for the reader's convenience. Here 

k-\-r-\-i 



k OO 



(/It) 



m— 1 



(2.41) 



r=0 



Tn=0 



r=0 



rl{k + r-iy. (fc + r)!(r-i)! 



By means of the above formula in the next theorem we derive an explicit expression for the state probabilities p^(i), 
k>0,t>0,a€ (0,1]. 

Theorem 2.4. The state probabilities p'^{t) = Pr{iV"(i) = A:|iV"(0) ^ i}, k > 0, t > 0, a e (0, 1], read 



A\ /A 



k OO k+r+i , ^ 

1 — n -m — n V / 



r— m— 
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E 

r=0 



k + 2r — i\ fk + 2r — i 
k + r 



where 



w'T((5 + r) 



^^Hr(/3r + 7) ^Hr(/3r + 7)r(<5)^ 

is i/ie Generalized Mittag-Lejjier function [12]. 

Proof. Recurring to Theorem 2.2 we can write for fc > 0, A 7^ /i, 



u;,7,/?,(5eC, 5R(/3) > 0, 



pm= / pi(s)Pr{ii;"(t)eds} 

Jo 



E 



A; 00 fc+r+'i 

E E 

r— m— 

1 



"Iml 



^ Vr!(A; + r - i)! (fc + r)!(r-i)! 



g-(A+M)Sgr-+m-lpj.|£,a(^) e ds} 



-(\+ti)s k+2r-i 



Applying the Laplace transform to both terms on the right-hand side we obtain 



A\ /A 



k 00 k-\-r-{-i 



E 



r— m— 

1 



^^\rl{k + r - iy. (fc + r)!(r-i)! 



A; CX3 /e+r+i 

E E 

1 



^Vr!(fc + r-i)! (fc + r)!(r-i)! 



-AV^-i 



00 

e-(>+/^)^s'-+™-i2"-ie'''^°ds 



z 



a-1 



\m\ (z" + A + ^iy+-^ 



(r + m — 1)! 



E 



A\ /A 



z-' + 



E E 



A Vt ^ — .^A k + 2r-iy. 



E 



r— m— 

k + 2r — i\ fk + 2r — i 



r — m I r + m 
r 



r 



k + r 



r + m 



(z" + A + /Lt)'-+™ 



(Z" + A + ^)fe+2r-i+l ■ 



(2.43) 



(2.44) 



(2.45) 



To invert equation (2.45) we use the Laplace transform (see formula (2.3.24) of Mathai and Haubold [13]) 
which immediately leads to (2.42). 

Remark 2.3. When a = l, formula (2.42) becomes the classical solution (2.41) because ^{w) — e^'/r{S) 



(2.46) 

□ 
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Remark 2.4. Result (2.42) is particularly interesting because its first addend contains the steady-state solution 



p^(i)=limKW= 1-- - , ^>0. (2.47) 



Furthermore, it is worth noticing that this geometric distribution coincides with that of the classical case a = I. The 
whole difference between the fractional and the non-fractional case lies in the transient regime. 

3 Path simulation and parameter estimation for the fractional linear 
birth-death process 

Wc now focus on a related point process which is relatively simpler to treat. Let J^{t), i > be a classical 
linear birth-death process with Afc > 0, /ifc > as its birth and death rates, respectively. Furthermore, define St, 
k gNU {0}, as the sojourn time of the process Af{t), t> in state k, i.e., given that the process is in state k, Sk is 
the time until the process leaves that state. It is well-known [18, Section 3.2, Chapter VI] that 

Pr{Sk>t} = exp[-iX + fi)kt], (3.1) 

and thus, 

Pr{Sk£dt}/dt^{\ + n)kexp[-{X + n)kt], t>0. (3.2) 

Analogous to the preceding section, in order to produce a fractional process related to the classical birth-death process 
it would be natural to substitute the unit-order time- derivative in the governing equations with a fractional derivative. 
This has been already carried out in Orsingher and Polito [14]. In the following, we exploit a subordination relation 
similar to that used in Section 2 in order to continue the analysis. In particular our aim is to develop methods 
suitable to simulation and parameter estimation that will be also applied to the fractional M/M/1 case in the last 
section. 

Recall thus that a fractional linear birth-death process A/'"(t), t > 0. a £ {0, 1] satisfies the subordination-relation 
[14] 

A/'"(i) = N[E"{t)], (3.3) 

where E°'{t) is the right- inverse process to an a-stable subordinator defined in the previous section. Using the above 
relation, we can easily calculate the distribution of the sojourn or holding times S*^, A: e N U {0} for the fractional 
linear birth-death process A/'"(i), t > 0, as follows. 

POO 

Pr{S'^ e dt)/dt = (A + ii)k / cxp [-(A + ^l)ks\ Pr{£;"(t) € ds} (3.4) 
= (A + y)kt°'-^Eo,,a [-(A + ^l)kt°'] , i > 0, fc > 1, < a < 1. 

Hence, the holding or sojourn time S*^ of the fractional linear birth-death process J\f°'{t) is Mittag-Leffler distributed. 
Another equivalent way to derive the event time distribution above is to replace the unit-order derivative in equation 
(3.4) of Taylor and Karlin [18, page 356] by the Caputo's fractional derivative operator D" used by [14]. That is, 

A" MS^ > - -(A + M)fcPr{5," > t}, (3.5) 
and solving the above equation, we obtain 

MSk >t} = [-(A + A*)fct"] , (3.6) 
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which gives equation (3.4). 

An interesting observation is that the birth and death sojourn times Bk and Dk-, respectively are no longer independent 
for < a < 1, i.e., 

P{Sk >t) = P(min{Bfe, > t) (3.7) 
^P{Bk>t,Dk>t) (3.8) 
^ P{Bk > t)P{Dk > t), (3.9) 

When the process is in state k, k G NU {0}, it transitions to the neighboring states k + 1 and k — 1 with probabilities 
A/(A + fj.) and At/(A + ^), respectively. Following Taylor and Karlin [18, page 358], a standard procedure to simulate 
trajectories of a fractional linear birth-death process is as follows: 

ALGORITHM : 

i) Fix the birth intensity A, the death intensity fi, and the initial population size A/'"(0) — m. 

ii) Simulate Sf = Sl^^T^ and U ^ U{0, 1). 



m 



liU < then Af"{si) = m + 1. Otherwise, A/'"(si) = m - 1. 



iv) Continuing in the same fashion and supposing that the current process state is Af"{sk-i) = k, generate 
d gi/'^r^^ and U = U{0, 1). lfU<j^ then A/'"(sfe) = k + 1. Otherwise, 7V"(sfc) = k - I. Repeat iv) 
until the desired population size is achieved or until extinct. 

Note that £k is negative exponentially distributed {£k = exp ((A + ^)fc)), and Tq, is a one-sided a-stable and 
independently distributed random variable. Below are some sample paths of the fractional linear birth-death process. 
The nonlinear trend, the longer holding times, and the slow or bursting behavior of the fractional linear birth-death 
process are apparent in Figure 1. 

We now provide point estimation algorithms for the parameters a, A, and /z. Assume that a sample trajectory of 
population size n corresponding to the n random inter-event times 5'^'s of the fractional linear birth-death process 
is observed, where there are ub births, no deaths, and n — hb + riB. Recall the structural representation of the 
Mittag-Leffler distributed random sojourn time = £^^"Tq,, where = exp {9k) is independent of a one-sided 
Q;+-stable distributed random variable Tq,, and 6 — X + Let = ln(5'£'). Then it is well-known that the mean 
and variance [see details in 7] of the log-transformed fc-th random sojourn time of the fractional linear birth-death 
process are 

- In (Ok) 

and 

respectively, where 7 ~ 0.5772156649 is the Euler-Mascheroni's constant. Following [6], the first two moments above 
therefore suggest that the simple linear regression model below can be fitted: 

S'^' = bo + blink + ek, fc = l,...,n, (3.12) 

where 

&o = ^^-7, 6i = -l/a, (3.13) 
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0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 



Figure 1: Sample paths of the fractional linear birth-death process: (top left) a = 1, A = 5,/i = 15; (top right) 
a = 0.8, A = 5,/i = 15; (bottom left) a = 0.8, A — 5,fi = 15; (bottom right) a — 0.8, A = /i = 5 with an initial 
population of m — 500. 



and Ek == (me = 0, 0"^ = o"^^, ) = In (^£^''°'Ta) +j,£= cxp(l). We point out that the error distribution depends 



only on a and is independent of the state k and 9, and this gives us a simple way of testing the rate fit as follows. 
Generate, say m samples (each of sample size n) from the error distribution using a. For a fixed significance level, 
test equality of two parent populations of the observed residuals and each of the m simulated errors using the 
two-sample Kolmogorov-Smirnov test, for instance. The proportion of the null acceptance out of to tests can then 
be used as a fit measure. 

By posing or , in formula (3.11) equal to its unbiased estimator ct^ = J2l=i ~ 2), we readily obtain the 

residual-based point estimators 



3 ^W-' 



-1/2 



exp 



(^-a (bo + 7) ) , 



(3.14) 



of the model parameters a and 9, correspondingly, where Sk = S"^ ~ ^ ^^"^ = ^0 + &i In k. Notice that the 
above estimators exploit the residuals to estimate a instead of the negative inverse of the least squares (LS) estimate 
of the slope 61. Furthermore, the least squares estimators of the slope and intercept are 



E;=i^/ (Inj-lnfc) 



60 = S?J -bi -Infc, 



(3.15) 



n 

where Infc = ^ \nj /n and 5^' ~ ^ S^/n. Hence, the closed-form point estimators of the intensities A and /i are 

^ of births 



A 



n 



(3.16) 
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and 

# of deaths A_ "^b^^^ , 



ji=- e = — ■e = e-x = — ■e={i — -]-e, (3.17) 

n n n \ n / 

respectively. Table 1 below shows some test results based on the percent bias 

[average estimate-parameter value] 

iUU X - 

parameter value 

and the coefficient of variation 

^ standard deviation of the estimates 
CV = 100 X 



average estimate 

using 1000 simulation runs. Note that we replaced the least squares estimator bo by the average of — (1/S) Infc 
to improve small sample performance. Apparently, the proposed point estimators, especially a performed relatively 
well even if the sample size is as small as 100. 



Table 1: Percent bias and dispersion of the proposed point estimators of a, A and fi. 



(a, A,/i) 


Estimator 


n = 
Bias 


10^ 

CV 


n = 
Bias 


103 
CV 


n = 
Bias 


10^ 
CV 




a 


0.989 


8.902 


0.075 


2.880 


0.010 


0.903 


(0.1,0.5,9) 


A 


22.389 


46.798 


2.920 


25.643 


0.394 


10.729 






24.938 


46.696 


3.568 


23.146 


0.363 


9.358 




a 


1.125 


8.919 


0.019 


2.667 


0.028 


0.909 


(0.25,1,6) 


A 


22.480 


45.615 


4.749 


21.624 


0.020 


9.309 






23.273 


45.997 


4.309 


20.478 


0.171 


9.467 




a 


0.709 


8.279 


0.140 


2.474 


0.016 


0.822 


(0.5,5,5) 


A 


20.249 


42.057 


3.343 


21.170 


0.545 


8.655 






20.992 


41.142 


3.767 


21.090 


0.621 


8.637 




a 


0.316 


7.213 


0.109 


2.272 


0.004 


0.679 


(0.75,7,1) 


A 


10.372 


38.858 


1.945 


17.192 


0.286 


6.991 






10.205 


42.200 


2.443 


19.205 


0.204 


7.472 




a 


0.924 


5.386 


0.077 


1.764 


0.008 


0.532 


(0.95,10,0.5) 


A 


9.544 


28.672 


1.539 


13.572 


0.272 


5.186 






8.875 


41.985 


1.559 


19.671 


0.468 


7.193 



We now provide formulas for the interval estimators of the model parameters. It is worth emphasizing that the 
explicit expressions of the estimators can be utilized to obtain resampling-based interval estimates especially for 
relatively small sample sizes. It is shown in [7] that 

^(a-a)^Nl0,^ ^j. (3.18) 

and a residual-based (1 — e)100% confidence interval for a directly follows as 

/S2 (32 - 20S2_34) ^^^^^ 
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where z^ji is the (1 — e/2)th quantile of the standard normal distribution, and < e < 1. We wiU now show the 
asymptotic normaUty of the estimators A and /I. 



Theorem 3.1. Letp = nB/nandp = X/9. Then 



as n oo where 



Vn{X-X) ^N{0, -p)+ p'al) 



(32 - 20q^ - a'^ 
40 



1 Ink 



(3.20) 



(3.21) 



Proof. It can be deduced from [6] and the asymptotic property of a Bernoulh/binomial sampled proportion that 



iV(0,S) 



(3.22) 



as n — > oo, where the variance-covariance matrix S is defined as 

Pil-P) 
al 

Invoking a standard result on asymptotic theory, the two-dimensional Central Limit Theorem implies that 

V^{h{e^) ~ h(0)) A N (o,h(0)Tsh(0)) , 



(3.23) 



(3.24) 



where On = {p, 0)"^ , h is a mapping from — > M, h(x) is continuous in a neighborhood of G M^, h{p, 9) = p ■ 9, 
and h{p,9) — {9,p) . This concludes the proof. □ 

Theorem 3.2. Let q = 1 — p = n]j/n, p = X/9, and q = Ijl/9. Then 



: _ ^) ^ TV (0, 9'p{l - pf + (1 - pfal) 



(3.25) 



as n — > cx) . 



Proof. The proof directly follows from the preceding theorem except that here we consider h(p, 9) = [1 — p) ■ 9 and 

Hp,e)^{~9,{i~p)f. □ 

We can now approximate the (1 — e)100% confidence interval for A and as 



A ± • d\, 



and 



respectively, where 



■pq+p^aj 



(3.26) 



(3.27) 
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and 

^, = ^^MIS^IM. (3.28) 

We now calculate the coverage probabilities using sample sizes n — lO'^, 10^, 10'* and 10'^ simulations to test our 
interval estimators of A and only. Note that the interval estimator for v has already been shown to perform well in 
past related studies [see 7, 6] . Table 2 below clearly suggests that the coverage probabilities of the interval estimators 
start to approach the true confidence level when n is at least 1000. If a narrower interval and a larger coverage 
are preferred then our simulations demonstrated that the previous point estimate replacement can be used. Note 
that this replacement and the above simple fit testing schemes can be directly applied to the fractional birth and 
fractional death processes in [6] to enhance performance of the point and interval estimators as well. 



Table 2: Coverage probabilities of the proposed interval estimators of X and fi using a 95% confidence level. 



(a, A, /i) 


Parameter 


n = 10^ 


n = 10^ 


n = 10"* 


(0.1,10,90) 


A 

fi 


0.879 
0.888 


0.936 
0.937 


0.954 
0.955 


(0.25,70,30) 


A 


0.892 
0.876 


0.925 
0.926 


0.955 
0.958 


(0.5,50,50) 


A 


0.896 
0.894 


0.934 
0.933 


0.946 
0.947 


(0.7,5,95) 


A 


0.864 


0.922 


0.947 


fi 


0.882 


0.924 


0.950 


(0.95,20,80) 


A 


0.900 
0.910 


0.948 
0.950 


0.951 
0.959 



Overall, the tables above provide additional merit to the proposed point and interval estimators of the model 
parameters. Aside from the computational simplicity of the proposed parameter estimation methods, the fit of the 
rates can also be checked easily. 

4 Trajectory generation and parameter estimation for the fractional 
simple linear birth-death or M/M/1 queue process 

Let N{t), t > he a classical simple birth-death process with A > 0, /i > as its constant birth and death rates, 
respectively. Furthermore, define Sk, fc e NU {0}, as the sojourn time of the process N{t), t > in state k, i.e. given 
that the process is in state fc, Sk is the time until the process leaves that state. Then from the preceding sections, 
it can easily be deduced that the holding/sojourn time S'^'s of the fractional simple birth-death or M/M/1 queue 
Na{t) are independently and identically (IID) Mittag-Leffler distributed, i.e., 

Pt{S^ edt}/dt^{X + n)t"-^E^^a.[-{X + n)t°'], i>0, ae(0,l]. (4.1) 

Note that everything here immediately follows from the previous section except that £k = exp (9) are IID. Assume 
that a sample trajectory of population size n corresponding to the n IID random inter-event times Sj^'s of the 
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fractional simple birth-death or M/M/1 process is observed, where there are ub births, njj deaths, and n = ub +fi_D- 
From [7] , a method-of- moments estimator for a is 

and 

^ = cxp I - a (ES*^ + 7) ) (4.3) 



is an estimator for 9. Recall that the asymptotic normality of a follows from the earlier result (3.18). Table 3 shows 
some test results based on the percent bias and the coefhcient of variation CV using 1000 simulations. Clearly, the 
proposed point estimators performed even better than the ones in the linear case. 



Table 3: Percent bias and dispersion of the proposed point estimators of a, A and ji. 



(a, A,/i) 


Estimator 


n = 
Bias 


102 
CV 


n = 
Bias 


10^ 
CV 


n = 
Bias 


lO'' 
CV 




a 


1.180 


9.077 


0.262 


2.924 


0.015 


0.864 


(0.1,0.5,9) 


A 


4.551 


45.894 


0.731 


15.929 


0.072 


4.804 




V- 


6.017 


25.210 


0.964 


8.649 


0.122 


2.728 




a 


0.973 


8.836 


0.044 


2.786 


0.022 


0.916 


(0.25,1,6) 


A 


7.344 


31.717 


0.270 


11.330 


0.054 


3.501 




V- 


5.322 


22.923 


0.337 


7.979 


0.021 


2.450 




a 


0.506 


8.043 


0.016 


2.534 


0.044 


0.818 


(0.5,5,5) 


A 


6.411 


24.339 


0.088 


8.609 


0.050 


2.640 




V- 


5.242 


28.822 


0.154 


8.869 


0.048 


2.682 




a 


0.793 


7.568 


0.054 


2.388 


0.027 


0.701 


(0.75,7,1) 


A 


3.933 


21.190 


0.452 


6.796 


0.204 


1.901 




V- 


6.330 


32.114 


0.497 


10.516 


0.098 


3.254 




a 


0.541 


5.587 


0.019 


1.795 


0.007 


0.558 


(0.95,10,0.5) 


A 


1.820 


13.947 


0.030 


4.508 


0.031 


1.540 




V- 


0.061 


43.363 


0.147 


14.001 


0.121 


4.539 



As in the preceding section, we provide formulas for the interval estimators of the model parameters. The explicit 
expressions of the estimators can be used to obtain resampling-based interval estimates especially for small sample 
sizes. A (1 — e)100% confidence interval for a is directly obtained from the previous section by simply replacing the 
point estimator of a. 

Theorem 4.1. Let p = ub /n and p = X/d . Then 



— >■ 00 where 



Vn(\-\] ^N{0, e'p{l -p)+ p^al) 



207r4(2 - a^) - i-n'^^a'' -f 20^2 - 32)(ln6l)2 - 720a3(ln 6')C(3) 



120^2 

and where (^(3) is the Riemann-zeta function evaluated at 3. 



(4.4) 



(4.5) 
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Proof. We omit the routine proof as it follows from the previous theorem. 



□ 



Theorem 4.2. Let q = 1 — p = tid/h, p = \/9, and q — n/O. Then 

^{r^-pi)^N (0, e^p{l - pf + (1 - pfal) (4.6) 

as n — > cxD . 

Proof. This directly follows from the preceding theorem. □ 

We now test our interval estimators by calculating the coverage probabilities using sample sizes n = 10^, 10^, 10^ and 
10^ simulations to test our interval estimators for A and Table 4 clearly suggests that the coverage probabilities 
of the interval estimators start to approach the true confidence level when n is at least 1000. 

Table 4: Coverage probabilities of the proposed interval estimators of A and fi using a 95% confidence level. 



(a, A,/i) 


Parameter 


n = 10^ 


n 10^ 


n = 10^ 


(0.1,10,90) 


A 


0.921 
0.932 


0.950 
0.959 


0.951 
0.948 


(0.25,70,30) 


A 


0.936 
0.927 


0.954 
0.942 


0.955 
0.958 


(0.5,50,50) 


A 


0.932 
0.925 


0.957 
0.948 


0.949 
0.953 


(0.7,5,95) 


A 


0.869 


0.944 


0.947 




0.931 


0.955 


0.950 


(0.95,20,80) 


A 


0.927 
0.917 


0.961 
0.947 


0.959 
0.947 



In general, the tables above indicate better performance of the proposed point and interval estimators than the 
procedures for the fractional linear birth-death process. 

5 Application 

We demonstrate our methods using the daily S&P index changes from January 1, 2011 until January 23, 2013. In 
particular, we apply the simple fractional birth-death or M/M/1 queue to model the number of positive- negative 
index changes. During the considered time period, there is a total of 516 events, of which 290 and 226 are positive 
and negative changes, respectively. Below is the trajectory of the number of the 516 positive-negative S&P index 
changes. 

Below are the point and 95% interval estimates for the S&P data, which seemingly indicates that the number of 
adjusted S&P index changes behaves like a standard simple birth-death process {a = 1). 

We examined the rate fit by simulating 1000 samples using the point estimates, and tested the equality of two parent 
populations using the two-sample Kolmogorov-Smirnov's test. Below are the summary statistics of the 1000 p-values. 

Note that the proportion of p-values larger than 0.05 is 99.9% which significantly indicates good model fit. The 
results above clearly indicated that the simple fractional birth-death model can also be used as a proof-of-concept or 
a smoothing-like tool for the standard birth-death process. 
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;ure 2: Path of the number of positive-negative S&P changes from 1/1/11 - 1/23/13. 



Table 5: Point and 95% interval estimates for the S&P data. 



Parameter 


Point estimate 


Interval estimate 


a 


0.995 


(0.950, 1.041) 


A 


0.053 


(0.042, 0.064) 




0.045 


(0.036, 0.054) 



Table 6: Summary statistics for the 1000 p-values. 



Min 


Qi 


Median 


Mean 


Q3 


Max 


0.033 


0.580 


0.786 


0.724 


0.942 


1.000 
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